{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d5371e04",
   "metadata": {},
   "source": [
    "# fenicsx的索引矩阵读写操作"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48d9c235",
   "metadata": {},
   "source": [
    "## 构建二维Poisson方程算例 \n",
    "使用P2元进行测试\n",
    "$$-\\Delta u = f\\ in\\ \\Omega\\\\\n",
    "u=g\\ on\\ \\partial\\Omega $$\n",
    "where $u=xe^y+ye^x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "ee8751ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "from dolfinx import mesh\n",
    "from mpi4py import MPI\n",
    "from dolfinx import fem, default_real_type\n",
    "import ufl\n",
    "import numpy as np\n",
    "from basix.ufl import element\n",
    "from dolfinx.fem.petsc import LinearProblem\n",
    "\n",
    "\n",
    "# 创建单位正方形网格，划分为nx×ny个四边形，再自动三角剖分\n",
    "nx, ny = 2, 2  # 4×4网格，h=0.25\n",
    "domain = mesh.create_unit_square(\n",
    "    MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle)\n",
    "p2_l = element(\"Lagrange\", domain.basix_cell(),\n",
    "               degree=2, dtype=default_real_type)\n",
    "V = fem.functionspace(domain, p2_l)  # 速度空间"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7031a86e",
   "metadata": {},
   "source": [
    "只要你的操作需要依赖网格实体之间的连接关系，就必须先调用`create_connectivity`。具体包括：\n",
    "* 访问单元格的顶点、边等组成元素\n",
    "* 定位边界或特定区域的网格实体\n",
    "* 施加边界条件（Dirichlet、Neumann 等）\n",
    "* 进行网格细化或自适应计算\n",
    "* 计算几何量（面积、体积、法向量等）\n",
    "* 在并行环境下处理分布式网格  \n",
    "\n",
    "`函数接口`\n",
    "```py\n",
    "(method) def create_connectivity(\n",
    "    d0: int,\n",
    "    d1: int\n",
    ") -> None\n",
    "Create connectivity between given pair of dimensions, d0 and d1.\n",
    "\n",
    "Args\n",
    "d0\n",
    "Dimension of entities one is mapping from\n",
    "\n",
    "d1\n",
    "Dimension of entities one is mapping to\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "c49a1145",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dolfinx.mesh.Mesh"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(domain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "187fe27b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 定义真解\n",
    "def u_exact(x):\n",
    "    return x[0] * np.exp(x[1]) + x[1] * np.exp(x[0])\n",
    "\n",
    "\n",
    "def f_exact(x):\n",
    "    return - (x[0] * np.exp(x[1]) + x[1] * np.exp(x[0]))\n",
    "\n",
    "\n",
    "def g_boundary(x):\n",
    "    return u_exact(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "c73fd3b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def search_boundary_edge_cell(msh: mesh.Mesh, counterclock=False) -> dict:\n",
    "    \"\"\"\n",
    "    识别并返回有限元网格中所有边界边及其对应的单元索引和坐标\n",
    "    \n",
    "    参数:\n",
    "        msh (mesh.Mesh): DOLFINx 网格对象\n",
    "        counterclock (bool): 是否按逆时针排序边界边，默认为 False\n",
    "    \n",
    "    返回:\n",
    "        dict: 字典，键为边界边的全局索引，值为包含单元索引和边坐标的字典\n",
    "    \"\"\"\n",
    "    # 获取网格拓扑维度信息\n",
    "    tdim = msh.topology.dim  # 网格维度 (2D=2)\n",
    "    fdim = tdim - 1          # 边界维度 (2D=1)\n",
    "    msh.topology.create_connectivity(fdim, tdim)\n",
    "    msh.topology.create_connectivity(fdim, 0)  # 创建边到顶点的连接关系\n",
    "\n",
    "    # 获取边到单元的连接关系\n",
    "    facet_to_cells = msh.topology.connectivity(fdim, tdim)\n",
    "    facet_to_vertices = msh.topology.connectivity(fdim, 0)  # 边到顶点的连接关系\n",
    "\n",
    "    # 获取网格几何信息（顶点坐标）\n",
    "    x = msh.geometry.x  # 所有顶点的坐标\n",
    "\n",
    "    # 构建结果字典\n",
    "    boundary_edges = {}\n",
    "    for facet_index in range(facet_to_cells.num_nodes):\n",
    "        cells = facet_to_cells.links(facet_index)\n",
    "        if cells.size == 1:  # 边界边\n",
    "            # 获取边的顶点索引\n",
    "            vertices = facet_to_vertices.links(facet_index)\n",
    "            # 获取顶点坐标\n",
    "            coords = x[vertices][:, :tdim]  # 提取前tdim个坐标分量（通常是2D）\n",
    "\n",
    "            boundary_edges[facet_index] = {\n",
    "                'cell_index': cells[0],  # 连接的单元索引\n",
    "                'coordinates': coords    # 边的顶点坐标（形状：(2, tdim)）\n",
    "            }\n",
    "    \n",
    "    # 如果需要逆时针排序，进行后处理\n",
    "    if counterclock and boundary_edges:\n",
    "        # 计算所有边界边的质心作为参考点\n",
    "        all_coords = np.vstack([edge['coordinates'] for edge in boundary_edges.values()])\n",
    "        centroid = np.mean(all_coords, axis=0)\n",
    "        \n",
    "        # 对每条边界边，检查其方向是否为逆时针，否则反转\n",
    "        for edge_data in boundary_edges.values():\n",
    "            p1, p2 = edge_data['coordinates']\n",
    "            # 计算向量从质心到p1和p2\n",
    "            v1 = p1 - centroid\n",
    "            v2 = p2 - centroid\n",
    "            # 计算叉积判断方向（z分量）\n",
    "            cross_product = v1[0] * v2[1] - v1[1] * v2[0]\n",
    "            # 如果叉积为负，说明是顺时针方向，需要反转\n",
    "            if cross_product < 0:\n",
    "                edge_data['coordinates'] = np.array([p2, p1])\n",
    "    \n",
    "    return boundary_edges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "d0b8dfab",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "网格中有8个单元格, nx = 2, ny = 2\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{0: {'cell_index': 0,\n",
       "  'coordinates': array([[0.5, 0. ],\n",
       "         [1. , 0. ]])},\n",
       " 3: {'cell_index': 2,\n",
       "  'coordinates': array([[0. , 0. ],\n",
       "         [0.5, 0. ]])},\n",
       " 4: {'cell_index': 0,\n",
       "  'coordinates': array([[1. , 0. ],\n",
       "         [1. , 0.5]])},\n",
       " 6: {'cell_index': 3,\n",
       "  'coordinates': array([[1. , 0.5],\n",
       "         [1. , 1. ]])},\n",
       " 11: {'cell_index': 4,\n",
       "  'coordinates': array([[0. , 0.5],\n",
       "         [0. , 0. ]])},\n",
       " 12: {'cell_index': 5,\n",
       "  'coordinates': array([[1. , 1. ],\n",
       "         [0.5, 1. ]])},\n",
       " 14: {'cell_index': 7,\n",
       "  'coordinates': array([[0. , 1. ],\n",
       "         [0. , 0.5]])},\n",
       " 15: {'cell_index': 7,\n",
       "  'coordinates': array([[0.5, 1. ],\n",
       "         [0. , 1. ]])}}"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 获取网格维度（2D问题为2）\n",
    "tdim = domain.topology.dim  # 2\n",
    "\n",
    "# 获取所有单元格的数量\n",
    "num_cells = domain.topology.index_map(tdim).size_local\n",
    "print(f\"网格中有{num_cells}个单元格, nx = {nx}, ny = {ny}\")\n",
    "\n",
    "# 读写网格的边界单元\n",
    "tdim = 2 # 拓扑维度 (2D网格=2)\n",
    "fdim = 1  # 边界维度(2D网格=1)\n",
    "boundary_facets = mesh.locate_entities_boundary(\n",
    "    domain, fdim,\n",
    "    lambda x: np.isclose(x[0], 0) | np.isclose(x[0], 1) |\n",
    "    np.isclose(x[1], 0) | np.isclose(x[1], 1)\n",
    ")  # 读取边界边索引\n",
    "\n",
    "search_boundary_edge_cell(domain, True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "9baa6b48",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 定义弱形式\n",
    "u = ufl.TrialFunction(V)\n",
    "v = ufl.TestFunction(V)\n",
    "f = fem.Function(V)\n",
    "f.interpolate(f_exact)  # 将右端项插值到函数空间\n",
    "a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx\n",
    "L = ufl.inner(f, v) * ufl.dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "4dc7c50b",
   "metadata": {},
   "outputs": [],
   "source": [
    "fdim = domain.topology.dim - 1  # 二维问题中边界是边，维度为 1\n",
    "boundary_facets = mesh.locate_entities_boundary(\n",
    "    domain, fdim,\n",
    "    lambda x: np.isclose(x[0], 0) | np.isclose(x[0], 1) |\n",
    "    np.isclose(x[1], 0) | np.isclose(x[1], 1)\n",
    ")\n",
    "boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)\n",
    "u_bc = fem.Function(V)\n",
    "u_bc.interpolate(g_boundary)\n",
    "bc = fem.dirichletbc(u_bc, boundary_dofs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "be82fe62",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nx = ny = 2, h = 0.5\n",
      "L2 norm = 0.0017322914731154812\n"
     ]
    }
   ],
   "source": [
    "# 求解\n",
    "problem = LinearProblem(a, L, bcs=[bc], petsc_options={\"ksp_type\": \"preonly\", \"pc_type\": \"lu\"})\n",
    "uh = problem.solve()\n",
    "# 真解\n",
    "u_true = fem.Function(V)\n",
    "u_true.interpolate(u_exact)\n",
    "\n",
    "# 打印误差\n",
    "error = u_true.x.array - uh.x.array\n",
    "# print(error)\n",
    "print(f\"nx = ny = {nx}, h = {1/nx}\")\n",
    "print(f\"L2 norm = {np.linalg.norm(error, ord=2)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a78a254b",
   "metadata": {},
   "source": [
    "h = 0.125 3.708008076609733e-05\n",
    "h = 0.0625 4.781744284825705e-6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8357449c",
   "metadata": {},
   "source": [
    "Poisson方程使用P2元 实现操作?  \n",
    "1. 获取单元格的所有节点索引\n",
    "2. 列出函数空间中自由度的坐标\n",
    "3. 给定参考坐标 映射为 物理坐标\n",
    "4. 给定物理坐标 计算对应基函数的值  \n",
    "3. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "33b0de27",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 操作1 获取单元格的所有节点索引（对P2元，每个单元有6个节点）\n",
    "cell_id = 1  #P2元共有 nx * ny * 2个单元格(三角形)\n",
    "dofmap = V.dofmap\n",
    "cell_nodes = dofmap.cell_dofs(cell_id)  # 返回该单元的6个节点的全局索引\n",
    "\n",
    "# 操作2 列出函数空间中自由度的坐标\n",
    "coordinates = V.tabulate_dof_coordinates()\n",
    "# shape = (2*nx+1)*(2*ny+1), 3 \n",
    "# 为兼容三维问题 这里第三个维度值为0\n",
    "# coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "7f882735",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "单元 0 的节点编号: [0 1 2 3 4 5]\n",
      "单元 1 的节点编号: [0 6 2 7 4 8]\n",
      "单元 2 的节点编号: [ 9  0  6  8 10 11]\n",
      "单元 3 的节点编号: [ 6  2 12 13 14  7]\n",
      "单元 4 的节点编号: [ 9 15  6 16 10 17]\n",
      "单元 5 的节点编号: [ 6 18 12 19 14 20]\n",
      "单元 6 的节点编号: [15  6 18 20 21 16]\n",
      "单元 7 的节点编号: [15 22 18 23 21 24]\n"
     ]
    }
   ],
   "source": [
    "for i in range (num_cells):\n",
    "    cell_node = dofmap.cell_dofs(i)  # 获取每个单元的节点索引\n",
    "    print(f\"单元 {i} 的节点编号: {cell_node}\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "0f496917",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.74326098, 1.61392435]])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 操作3\n",
    "# 给定参考坐标 映射为 物理坐标\n",
    "def compute_physical_point(p0, p1, p2, X):\n",
    "    \"\"\"\n",
    "    Map coordinates `X` in reference element to triangle defined by `p0`, `p1` and `p2`\n",
    "    \"\"\"\n",
    "    el = element(\"Lagrange\", \"triangle\", 1)\n",
    "    basis_values = el.tabulate(0, X)\n",
    "    x = basis_values[0] @ np.vstack([p0, p1, p2])\n",
    "    return x\n",
    "\n",
    "\n",
    "p0 = np.array([2.0, 1.8])\n",
    "p1 = np.array([1.0, 1.2])\n",
    "p2 = np.array([1.3, 1.0])\n",
    "reference_points = np.array([[0.19773226, 0.08429537]])  # 这里传入二维数组\n",
    "x = compute_physical_point(p0, p1, p2, reference_points)\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d0ac07b",
   "metadata": {},
   "source": [
    "## 操作4 \n",
    "给定物理坐标 计算对应基函数的值  \n",
    "数学过程是\n",
    "$$\n",
    "\\tilde{u}(x)=\\sum^n_{i=1}{N_i(x)\\cdot u_i}\n",
    "$$\n",
    "其中$N_i(x)$为节点对应的基函数, $u_i$为节点$i$处的物理量(有限元解$Ax=b$算出的结果), $n=6$(P2元)    \n",
    "fenicsx节点编号为  \n",
    "2  \n",
    "4 3  \n",
    "0 5 1  \n",
    "`tabulate`本质是一个语法糖 负责计算标准基函数下 坐标对应的基函数值"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c184531",
   "metadata": {},
   "source": [
    "![p2-basis-function.png](../assets/p2-basis-function.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "5f582b73",
   "metadata": {},
   "outputs": [],
   "source": [
    "basis_functions = [\n",
    "    lambda xi, eta: 2*xi**2 + 2*eta**2 + 4*xi*eta - 3*xi - 3*eta + 1,  # N1\n",
    "    lambda xi, eta: 2*xi**2 - xi,                                      # N2\n",
    "    lambda xi, eta: 2*eta**2 - eta,                                      # N3\n",
    "    lambda xi, eta: -4*xi**2 - 4*xi*eta + 4*xi,                        # N4\n",
    "    lambda xi, eta: 4*xi*eta,                                           # N5\n",
    "    lambda xi, eta: -4*eta**2 - 4*xi*eta + 4*eta                       # N6\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "189d6816",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_basis_values(xi, eta):\n",
    "    \"\"\"\n",
    "    Compute basis function values at given reference coordinates (xi, eta)\n",
    "    \"\"\"\n",
    "    return np.array([N(xi, eta) for N in basis_functions])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0faaa61b",
   "metadata": {},
   "source": [
    "## dolfinx.geometry\n",
    "Methods for geometric searches and operations.  \n",
    "本节的碰撞不是物理中的动力学碰撞，而是判断点、线段、面或体等几何实体是否与网格单元在空间中相交或者包含。  \n",
    "`bb_tree`将网格实体（如单元）用最小包围盒包裹，加速碰撞检测。例如，当判断点是否与某个单元 “碰撞” 时，先检查点是否在包围盒内，再进行精确几何计算，减少计算量。   \n",
    "`compute_collisions_points`：该函数利用包围盒树快速筛选与点碰撞的包围盒叶子节点，返回每个点对应的候选包围盒索引（用于后续精确检测的预处理）。\n",
    "`compute_colliding_cells`：该函数用于根据候选单元格列表和给定点集，精确计算与点发生几何碰撞（点在单元内或边界上）的网格单元索引。  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "77905bcd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.464\n",
      "物理坐标 [1 1] → 参考坐标 (ξ,η) = (0.2113, 0.5774)\n"
     ]
    }
   ],
   "source": [
    "def jacobian(p0, p1, p2):\n",
    "    \"\"\"计算雅可比矩阵的行列式（三角形面积的2倍）\"\"\"\n",
    "    x1, y1 = p0\n",
    "    x2, y2 = p1\n",
    "    x3, y3 = p2\n",
    "    return (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)\n",
    "\n",
    "\n",
    "def physical_to_reference(p0, p1, p2, physical_coords):\n",
    "    \"\"\"将物理坐标映射到参考坐标系 (ξ,η) ∈ [0,1]x[0,1] 且 ξ+η ≤ 1\"\"\"\n",
    "    x, y = physical_coords\n",
    "    x1, y1 = p0\n",
    "    x2, y2 = p1\n",
    "    x3, y3 = p2\n",
    "\n",
    "    # 计算三个子三角形的面积\n",
    "    L1 = jacobian(physical_coords, p1, p2)\n",
    "    L2 = jacobian(physical_coords, p2, p0)\n",
    "    L3 = jacobian(physical_coords, p0, p1)\n",
    "    J = jacobian(p0, p1, p2)\n",
    "\n",
    "    # 面积坐标转换为参考坐标\n",
    "    xi = L2 / J\n",
    "    eta = L3 / J\n",
    "\n",
    "    return xi, eta\n",
    "\n",
    "\n",
    "# 测试\n",
    "p0 = np.array([0, 0])\n",
    "p1 = np.array([2, 0])\n",
    "p2 = np.array([1, 1.732])\n",
    "print(jacobian(p0, p1, p2))\n",
    "physical_coords = np.array([1, 1])\n",
    "xi, eta = physical_to_reference(p0, p1, p2, physical_coords)\n",
    "print(f\"物理坐标 {physical_coords} → 参考坐标 (ξ,η) = ({xi:.4f}, {eta:.4f})\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "ee644264",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<AdjacencyList> with 1 nodes\n",
      "  0: [4 ]\n",
      "\n",
      "[0.22 0.24] FEM解: 0.5767235743307436\n",
      "[0.22 0.24] 真解: 0.5787332284116804\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "dolfinx.mesh.Mesh"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from dolfinx import geometry\n",
    "bb_tree = geometry.bb_tree(domain, 2)  # 创建边界包围盒树\n",
    "\n",
    "# 定义目标物理点\n",
    "target_point = np.array([0.22, 0.24, 0.0])  # 3D点 (z=0 for 2D)\n",
    "# 查找包含点的单元格\n",
    "candidates = geometry.compute_collisions_points(bb_tree, target_point)\n",
    "target_cell = geometry.compute_colliding_cells(\n",
    "    domain, candidates, target_point)\n",
    "print(target_cell)\n",
    "target_cell_id = target_cell.array[0]\n",
    "cell_nodes_idx = dofmap.cell_dofs(target_cell_id)  # 返回该单元的6个节点的全局索引\n",
    "\n",
    "# 获取cell_nodes_idx对应的有限元结果\n",
    "uh_values = uh.x.array[cell_nodes_idx]  # 获取该单元格的所有节点的有限元结果\n",
    "\n",
    "# V.tabulate_dof_coordinates().shape  # 返回空间的自由度坐标 shape = (2*nx+1)*(2*ny+1), 3\n",
    "# 获取单元格的所有节点坐标\n",
    "coordinates = V.tabulate_dof_coordinates()[cell_nodes_idx]  # 获取该单元的所有节点的坐标\n",
    "\n",
    "p0 = coordinates[0][:-1]  # 第一个节点\n",
    "p1 = coordinates[1][:-1]  # 第二个节点\n",
    "p2 = coordinates[2][:-1]  # 第三个节点\n",
    "\n",
    "reference_point = physical_to_reference(p0, p1, p2, target_point[:2])\n",
    "basis_ref = p2_l.tabulate(0, np.array([reference_point])) \n",
    "# 参考域与物理域中对应节点的基函数值完全相同\n",
    "print(f\"{target_point[:2]} FEM解: {np.dot(basis_ref[0], uh_values)[0]}\")\n",
    "print(f\"{target_point[:2]} 真解: {u_exact(target_point[:2])}\")\n",
    "type(domain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad76e077",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "79875424",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.  , 0.  ],\n",
       "       [0.  , 0.5 ],\n",
       "       [0.5 , 0.5 ],\n",
       "       [0.25, 0.5 ],\n",
       "       [0.25, 0.25],\n",
       "       [0.  , 0.25]])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coordinates[:, :-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98959b10",
   "metadata": {},
   "source": [
    "# 矢量元的补0现象\n",
    "标量P1元 传入一个位于标准三角形内的坐标点 返回长度为3的基函数值列表    \n",
    "2D P1元(x方向 和 y方向) 传入一个位于标准三角形内的坐标点 返回长度为2*6的基函数值列表(实际尺寸为1\\*1\\*2\\*6)    \n",
    "存储值时，tabulate函数会将数组进行扩展，数组长度为2(P1元的维度)\\*3(P1元有三个基函数)=6  \n",
    "第一个列表只存储x方向的值 y方向的值计为0；\n",
    "第二个列表只存储y方向的值 x方向的值直接计为0；  \n",
    "可以观察下面的程序\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "72ce20df",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p1 矢量 = [[-2.77555756e-17  0.00000000e+00  5.00000000e-01  0.00000000e+00\n",
      "   5.00000000e-01  0.00000000e+00]\n",
      " [ 0.00000000e+00 -2.77555756e-17  0.00000000e+00  5.00000000e-01\n",
      "   0.00000000e+00  5.00000000e-01]]\n",
      "p1 标量 = [[[-2.77555756e-17  5.00000000e-01  5.00000000e-01]]]\n"
     ]
    }
   ],
   "source": [
    "# 矢量元 tabulate函数\n",
    "# coord = np.random.rand(2) / 2\n",
    "coord = [0.5, 0.5]\n",
    "solid_msh = mesh.create_rectangle(\n",
    "    MPI.COMM_WORLD, [[-1, -1], [1, 1]], n=[4, 4], cell_type=mesh.CellType.triangle)\n",
    "lambda_p1 = element(\n",
    "    \"Lagrange\", solid_msh.basix_cell(), degree=1, shape=(solid_msh.geometry.dim,), dtype=default_real_type\n",
    ")\n",
    "lambda_p1_1 = element(\"Lagrange\", solid_msh.basix_cell(),\n",
    "               degree=1, dtype=default_real_type)  # 标量p1元\n",
    "basis_values = lambda_p1.tabulate(0, np.array([coord]))\n",
    "bb = lambda_p1_1.tabulate(0, np.array([coord]))\n",
    "print(f\"p1 矢量 = {basis_values[0][0]}\")\n",
    "print(f\"p1 标量 = {bb}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "47f4ded9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[-2.77555756e-17,  5.00000000e-01,  5.00000000e-01]],\n",
       "\n",
       "       [[-1.00000000e+00,  1.00000000e+00,  0.00000000e+00]],\n",
       "\n",
       "       [[-1.00000000e+00,  0.00000000e+00,  1.00000000e+00]]])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rs = lambda_p1_1.tabulate(1, np.array([coord]))\n",
    "rs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "6c21db88",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 3 5 7]\n"
     ]
    }
   ],
   "source": [
    " # 选取奇数索引\n",
    "a = np.array([1, 2, 3, 4, 5, 6, 7])\n",
    "b = a[::2]\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "450ea264",
   "metadata": {},
   "source": [
    "## basix.quadrature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "96c16655",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "积分点坐标: [[0.16666667 0.16666667]\n",
      " [0.16666667 0.66666667]\n",
      " [0.66666667 0.16666667]]\n",
      "积分权重: [0.16666667 0.16666667 0.16666667]\n"
     ]
    }
   ],
   "source": [
    "import basix\n",
    "from basix import CellType, QuadratureType\n",
    "\n",
    "# 获取积分点和权重（直接解包元组）\n",
    "q_points, q_weights = basix.make_quadrature(\n",
    "    CellType.triangle,  # 位置参数 cell\n",
    "    degree=2,           # 关键字参数 degree\n",
    "    rule=QuadratureType.default  # 枚举成员小写\n",
    ")\n",
    "\n",
    "    \n",
    "print(\"积分点坐标:\", q_points)\n",
    "print(\"积分权重:\", q_weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "fea43a59",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "积分结果: 0.25\n",
      "估计误差: 5.527085942945804e-15\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import dblquad\n",
    "\n",
    "# 定义被积函数\n",
    "\n",
    "\n",
    "def f(x, y):\n",
    "    return x**2 + y # 示例函数\n",
    "\n",
    "\n",
    "# 计算在标准三角形上的积分\n",
    "# x ∈ [0,1], y ∈ [0, 1-x]\n",
    "result, error = dblquad(\n",
    "    lambda y, x: f(x, y),  # 注意：dblquad 参数顺序是 (y, x)\n",
    "    0, 1,                  # x 下限\n",
    "    lambda x: 0,           # y 下限函数\n",
    "    lambda x: 1 - x        # y 上限函数\n",
    ")\n",
    "\n",
    "print(f\"积分结果: {result}\")\n",
    "print(f\"估计误差: {error}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26596e13",
   "metadata": {},
   "source": [
    "## gauss积分测试"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "fec369d9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "积分结果: 4.38905609893065\n"
     ]
    },
    {
     "ename": "AttributeError",
     "evalue": "module 'tools' has no attribute 'jacobian'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[42], line 30\u001b[0m\n\u001b[1;32m     23\u001b[0m vortex \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39marray([\n\u001b[1;32m     24\u001b[0m     [\u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m0\u001b[39m],\n\u001b[1;32m     25\u001b[0m     [\u001b[38;5;241m2\u001b[39m, \u001b[38;5;241m0\u001b[39m],\n\u001b[1;32m     26\u001b[0m     [\u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m2\u001b[39m],\n\u001b[1;32m     27\u001b[0m ])\n\u001b[1;32m     28\u001b[0m phy_coords \u001b[38;5;241m=\u001b[39m tools\u001b[38;5;241m.\u001b[39mreference_to_physical(\n\u001b[1;32m     29\u001b[0m     lambda_p1_1, \u001b[38;5;241m*\u001b[39mvortex, tools\u001b[38;5;241m.\u001b[39mGaussIntegrate()\u001b[38;5;241m.\u001b[39mpoints)\n\u001b[0;32m---> 30\u001b[0m det_J \u001b[38;5;241m=\u001b[39m \u001b[43mtools\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mjacobian\u001b[49m(\u001b[38;5;241m*\u001b[39mvortex)\n\u001b[1;32m     31\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdet_J = \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mdet_J\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m     32\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mphy_coords = \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mphy_coords\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n",
      "\u001b[0;31mAttributeError\u001b[0m: module 'tools' has no attribute 'jacobian'"
     ]
    }
   ],
   "source": [
    "from scipy import integrate\n",
    "import tools\n",
    "import numpy as np\n",
    "from basix.ufl import element\n",
    "def f1(x):\n",
    "    return x * np.exp(2*x) - x\n",
    "\n",
    "\n",
    "def f2(x):\n",
    "    return x * np.exp(-2*x+4) - x\n",
    "\n",
    "def ff(x, y):\n",
    "    return x * np.exp(y)\n",
    "\n",
    "bd1 = [0, 1]\n",
    "bd2 = [1, 2]\n",
    "# 计算数值积分\n",
    "r1 = integrate.quad(f1, 0, 1)\n",
    "r2 = integrate.quad(f2, 1, 2)\n",
    "print(f\"积分结果: {r1[0] + r2[0]}\")\n",
    "\n",
    "# tools gs \n",
    "vortex = np.array([\n",
    "    [0, 0],\n",
    "    [2, 0],\n",
    "    [1, 2],\n",
    "])\n",
    "phy_coords = tools.reference_to_physical(\n",
    "    lambda_p1_1, *vortex, tools.GaussIntegrate().points)\n",
    "det_J = tools.jacobian(*vortex)\n",
    "print(f\"det_J = {det_J}\")\n",
    "print(f\"phy_coords = {phy_coords}\")\n",
    "# gs\n",
    "ff(phy_coords[:, 0], phy_coords[:, 1])\n",
    "gs_res = det_J * \\\n",
    "    np.dot(ff(phy_coords[:, 0], phy_coords[:, 1]),\n",
    "           tools.GaussIntegrate().weights)\n",
    "print(f\"gs_res = {gs_res}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f8b3ce1",
   "metadata": {},
   "source": [
    "## P2元变量存储顺序\n",
    "是交错排列还是先排x分量后排列y分量  \n",
    "验证不难发现是交错排列的..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d42c4da9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from mpi4py import MPI\n",
    "from basix.ufl import element\n",
    "from dolfinx import mesh, default_real_type\n",
    "import ufl\n",
    "\n",
    "def true_velocity(x):\n",
    "    return [x[0], x[1]+10]\n",
    "\n",
    "n_list = [1, 1]\n",
    "msh = mesh.create_rectangle(\n",
    "    MPI.COMM_WORLD, [np.array([0, 0]), np.array(\n",
    "        [1, 1])], n_list, mesh.CellType.triangle\n",
    ")\n",
    "x = ufl.SpatialCoordinate(msh)\n",
    "P2 = element(\n",
    "    \"Lagrange\", msh.basix_cell(), degree=2, shape=(msh.geometry.dim,), dtype=default_real_type\n",
    ")\n",
    "V = fem.functionspace(msh, P2)\n",
    "u_true = fem.Function(V)\n",
    "u_true.interpolate(true_velocity)\n",
    "u_true.x.array\n",
    "# array([0., 10.,  1., 10.,  1., 11.,  1., 10.5,  0.5, 10.5,  0.5,\n",
    "    #    10.,  0., 11.,  0.5, 11.,  0., 10.5])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94f01c2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1, 2, 3, 4])\n",
    "a[-2:]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "fenicsx-env",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
